{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pseudo-compressible Navier-Stokes element symbolic formulation and implementation\n",
    "**Rubén Zorrilla, Riccardo Rossi**\n",
    "UPC BarcelonaTech, CIMNE\n",
    "\n",
    "The weakly-compressible Navier-Stokes problem is conformed by the linear momentum conservation equation\n",
    "\\begin{equation*}\n",
    "\\frac{\\partial\\rho\\mathbf{u}}{\\partial t} + \\rho\\mathbf{a}\\cdot\\nabla\\mathbf{u} - \\nabla\\cdot\\mathbf{\\boldsymbol{\\sigma}} + \\frac{\\mathbf{u}}{\\kappa} = \\rho\\mathbf{f}\n",
    "\\end{equation*}\n",
    "together with the mass conservation equation\n",
    "\\begin{equation*}\n",
    "\\frac{D\\rho}{Dt} + \\rho\\nabla\\cdot\\mathbf{u} = 0\n",
    "\\end{equation*}\n",
    "with $\\mathbf{u}$ denoting the velocity field, $\\rho$ the density and $\\mathbf{f}$ the volumetric force.\n",
    "For the sake of generality a Darcy term governed by the $\\kappa$ the permeability is also included in the formulation.\n",
    "\n",
    "Developing the material time derivative of the mass conservation equation one obtains the equation\n",
    "\\begin{equation*}\n",
    "\\frac{\\partial\\rho}{\\partial t} + \\mathbf{u}\\cdot\\nabla\\rho + \\rho\\nabla\\cdot\\mathbf{u} = 0\\; ,\n",
    "\\end{equation*}\n",
    "which intentionally keeps the $\\nabla \\rho$ term as it might be non-negligible in the weakly-compressible scenario.\n",
    "\n",
    "The density local time derivative can be rewritten in terms of the pressure time derivative as\n",
    "\\begin{equation*}\n",
    "\\frac{\\partial p}{\\partial t} = \\frac{\\partial p}{\\partial\\rho}\\frac{\\partial\\rho}{\\partial t} \\longrightarrow \\frac{\\partial\\rho}{\\partial t} = \\left(\\frac{\\partial p}{\\partial\\rho}\\right)^{-1}\\frac{\\partial p}{\\partial t}\\; .\n",
    "\\end{equation*}\n",
    "\n",
    "Introducing the state equation\n",
    "\\begin{equation*}\n",
    "\\frac{\\partial p}{\\partial\\rho} = c^{2}\n",
    "\\end{equation*}\n",
    "into the mass conservation equation yields\n",
    "\\begin{equation*}\n",
    "\\frac{1}{c^{2}}\\frac{\\partial p}{\\partial t} + \\mathbf{u}\\cdot\\nabla\\rho + \\rho\\nabla\\cdot\\mathbf{u} = 0\n",
    "\\end{equation*}\n",
    "where $c$ is the sound velocity of the considered fluid.\n",
    "\n",
    "Once arrived to this point, one can divide the previous equation by $\\rho$ for the sake of having a better conditioned matrix (specially in those cases in where $\\rho\\gg1$).\n",
    "This has been implemented as an user-defined option in the symbolic element generation script. Note that this will also affect the stabilization terms.\n",
    "\n",
    "The density is considered to be time independent in the momentum conservation equation, implying that\n",
    "\\begin{equation*}\n",
    "\\rho\\frac{\\partial\\mathbf{u}}{\\partial t} + \\rho\\mathbf{a}\\cdot\\nabla\\mathbf{u} - \\nabla\\cdot\\mathbf{\\boldsymbol{\\sigma}} + \\frac{\\mathbf{u}}{\\kappa} = \\rho\\mathbf{f}\\; .\n",
    "\\end{equation*}\n",
    "\n",
    "In addition, the stress tensor $\\boldsymbol{\\sigma}$ can be split in an isochoric and deviatoric (viscous) contribution as \n",
    "\\begin{equation*}\n",
    "\\boldsymbol{\\sigma} = -p\\mathbf{I} + \\rho(\\nabla\\cdot\\mathbf{u})\\mathbf{I} + \\boldsymbol{\\tau}\\; .\n",
    "\\end{equation*}\n",
    "Taking into account that no volumetric deformation is considered in the momentum conservation equation, $\\boldsymbol{\\sigma}$ can be simplified to\n",
    "\\begin{equation*}\n",
    "\\boldsymbol{\\sigma}=-p\\mathbf{I}+\\boldsymbol{\\tau}\\; .\n",
    "\\end{equation*}\n",
    "where $p$ is the pressure and $\\boldsymbol{\\tau}$ the shear stress tensor, which is obtained from the corresponding constitutive model as\n",
    "\\begin{equation*}\n",
    "\\boldsymbol{\\tau} = \\mathbb{\\mathbb{C}}\\nabla^{s}\\mathbf{u}\\; ,\n",
    "\\end{equation*}\n",
    "being $\\nabla^{s}$ the symmetric gradient operator and $\\mathbb{C}$ the 4th order constitutive tensor.\n",
    "\n",
    "Substituting the stress tensor above, the linear momentum equation reads\n",
    "\\begin{equation*}\n",
    "\\rho\\frac{\\partial\\mathbf{u}}{\\partial t} + \\rho\\mathbf{a}\\cdot\\nabla\\mathbf{u} - \\nabla\\cdot\\left(-p\\mathbf{I}+\\boldsymbol{\\tau}\\right) + \\frac{\\mathbf{u}}{\\kappa} = \\rho\\mathbf{f}\n",
    "\\end{equation*}\n",
    "\n",
    "Finally, the linear momentum and mass conservation residuals can be defined as\n",
    "\\begin{align*}\n",
    "\\mathbf{R}^{M}(\\mathbf{u},p) = & \\rho\\mathbf{f} - \\rho\\frac{\\partial\\mathbf{u}}{\\partial t} - \\rho\\mathbf{a}\\cdot\\nabla\\mathbf{u} + \\nabla\\cdot\\left(-p\\mathbf{I}+\\boldsymbol{\\tau}\\right) - \\frac{\\mathbf{u}}{\\kappa} \\\\\n",
    "\\mathbf{R}^{C}(\\mathbf{u},p) = & -\\frac{1}{c^{2}}\\frac{\\partial p}{\\partial t} - \\mathbf{u}\\cdot\\nabla\\rho - \\rho\\nabla\\cdot\\mathbf{u}\n",
    "\\end{align*}\n",
    "\n",
    "## Galerkin weak form and ASGS stabilization\n",
    "\n",
    "Considering the velocity test function $\\mathbf{w}$ and pressure test function $q$, the Galerkin weak form of the problem can be obtained as\n",
    "\n",
    "\\begin{equation*}\n",
    "\\left(\\mathbf{w},\\mathbf{R}^{M}(\\mathbf{u},p)\\right)_{\\Omega}+\\left(q,\\mathbf{R}^{C}(\\mathbf{u},p)\\right)_{\\Omega} = 0\\; ,\n",
    "\\end{equation*}\n",
    "that is\n",
    "\\begin{equation*}\n",
    "\\left(\\mathbf{w}, \\rho\\mathbf{f} - \\rho\\frac{\\partial\\mathbf{u}}{\\partial t} - \\rho\\mathbf{a}\\cdot\\nabla\\mathbf{u} + \\nabla\\cdot\\left(-p\\mathbf{I}+\\boldsymbol{\\tau}\\right) - \\frac{\\mathbf{u}}{\\kappa}\\right)_{\\Omega} + \\left(q, -\\frac{1}{c^{2}}\\frac{\\partial p}{\\partial t} - \\mathbf{u}\\cdot\\nabla\\rho - \\rho\\nabla\\cdot\\mathbf{u}\\right)_{\\Omega} = 0\\; .\n",
    "\\end{equation*}\n",
    "\n",
    "On top of that, the following multiscale decomposition is considered\n",
    "\n",
    "\\begin{align*}\n",
    "\\mathbf{u} & = \\mathbf{u_{h}}+\\mathbf{u_{s}}\\\\\n",
    "p & = p_{h}+p_{s}\n",
    "\\end{align*}\n",
    "\n",
    "where $\\mathbf{u_{\\mathbf{h}}}$ and $p_{h}$ belong to the finite element space while $\\mathbf{u_{s}}$ and $p_{s}$ are the so called unresolvable subscales. Substituting the proposed decomposition in the Galerkin residual weak form above yields\n",
    "\n",
    "\\begin{equation*}\n",
    "\\left(\\mathbf{w},\\rho\\mathbf{f} - \\rho\\frac{\\partial(\\mathbf{u}_{h}+\\mathbf{u}_{s})}{\\partial t} - \\rho\\mathbf{a}\\cdot\\nabla(\\mathbf{u}_{h}+\\mathbf{u}_{s}) + \\nabla\\cdot\\left(-(p_{h}+p_{s})\\mathbf{I}+\\mathbb{\\mathbb{C}}\\nabla^s(\\mathbf{u}_{h}+\\mathbf{u}_{s})\\right) - \\frac{\\left(\\mathbf{u}_{h}+\\mathbf{u}_{s}\\right)}{\\kappa} \\right)_{\\Omega} + \\left(q,-\\frac{1}{c^{2}}\\frac{\\partial(p_{h}+p_{s})}{\\partial t} - \\mathbf{a}\\cdot\\nabla\\rho - \\rho\\nabla\\cdot(\\mathbf{u}_{h}+\\mathbf{u}_{s})\\right)_{\\Omega} = 0\n",
    "\\end{equation*}\n",
    "\n",
    "At this point it is important to clearly state the assumptions concerning the subscales that are to be taken from now on\n",
    "- Quasi-static subscales: $\\frac{\\partial\\mathbf{u_{s}}}{\\partial t}=\\frac{\\partial p_{s}}{\\partial t}\\approx0$\n",
    "- $\\mathbf{u_{s}}\\approx p_{s}\\approx0$ on $\\Gamma$\n",
    "- Linear order finite elements implying that higher order derivatives are equal to 0\n",
    "\n",
    "Developing terms in the decomposed Galerkin residual weak form considering the quasi-static nature of the subscales one obtains\n",
    "\\begin{align*}\n",
    "\\left(\\mathbf{w},\\rho\\mathbf{f}\\right)_{\\Omega} &- \\left(\\mathbf{w},\\rho\\frac{\\partial\\mathbf{u}_{h}}{\\partial t}\\right)_{\\Omega} - \\left(\\mathbf{w},\\rho\\mathbf{a}\\cdot\\nabla\\mathbf{u}_{h}\\right)_{\\Omega} - \\left(\\mathbf{w},\\rho\\mathbf{a}\\cdot\\nabla\\mathbf{u}_{s}\\right)_{\\Omega} - \\left(\\mathbf{w},\\nabla p_{h}\\right)_{\\Omega} - \\left(\\mathbf{w},\\nabla p_{s}\\right)_{\\Omega} + \\left(\\mathbf{w},\\nabla\\cdot\\mathbb{\\mathbb{C}}\\nabla^s\\mathbf{u}_{h}\\right)_{\\Omega} + \\left(\\mathbf{w},\\nabla\\cdot\\mathbb{\\mathbb{C}}\\nabla^s\\mathbf{u}_{s}\\right)_{\\Omega} \\\\ \n",
    "&- \\left(\\mathbf{w},\\frac{\\mathbf{u}_{h}}{\\kappa}\\right)_{\\Omega} - \\left(\\mathbf{w},\\frac{\\mathbf{u}_{s}}{\\kappa}\\right)_{\\Omega} - \\left(q,\\frac{1}{c^{2}}\\frac{\\partial p_{h}}{\\partial t}\\right)_{\\Omega} - \\left(q,\\mathbf{a}\\cdot\\nabla\\rho\\right)_{\\Omega} - \\left(q,\\rho\\nabla\\cdot\\mathbf{u}_{h}\\right)_{\\Omega} - \\left(q,\\rho\\nabla\\cdot\\mathbf{u}_{s}\\right)_{\\Omega} = 0 \\; .\n",
    "\\end{align*}\n",
    "\n",
    "Integrating by parts the shear stress terms\n",
    "\\begin{align*}\n",
    "\\left(\\mathbf{w},\\nabla\\cdot\\mathbb{\\mathbb{C}}\\nabla^s\\mathbf{u}_{h}\\right)_{\\Omega} & = -\\left(\\nabla \\mathbf{w},\\mathbb{\\mathbb{C}}\\nabla^s\\mathbf{u}_{h}\\right)_{\\Omega} + \\langle\\mathbf{w},\\left(\\mathbb{C}\\nabla^s\\mathbf{u}_{h}\\right)\\cdot\\mathbf{n}\\rangle_{\\Gamma}\\; ,\\\\\n",
    "\\left(\\mathbf{w},\\nabla\\cdot\\mathbb{\\mathbb{C}}\\nabla^s\\mathbf{u}_{s}\\right)_{\\Omega} & = -\\left(\\nabla \\mathbf{w},\\mathbb{\\mathbb{C}}\\nabla^s\\mathbf{u}_{s}\\right)_{\\Omega} + \\langle\\mathbf{w},\\left(\\mathbb{C}\\nabla^s\\mathbf{u}_{s}\\right)\\cdot\\mathbf{n}\\rangle_{\\Gamma} \\approx - \\left(\\nabla\\mathbf{w},\\mathbb{\\mathbb{C}}\\nabla^s\\mathbf{u}_{s}\\right)_{\\Omega}\\; ,\n",
    "\\end{align*}\n",
    "and doing so for the pressure ones\n",
    "\\begin{align*}\n",
    "\\left(\\mathbf{w},\\nabla p_{h}\\right)_{\\Omega} & = -\\left(\\nabla\\cdot\\mathbf{w},p_{h}\\right)_{\\Omega} + \\langle\\mathbf{w}\\cdot\\mathbf{n}, p_{h}\\rangle_{\\Gamma}\\\\\n",
    "\\left(\\mathbf{w},\\nabla p_{s}\\right)_{\\Omega} & = -\\left(\\nabla\\cdot\\mathbf{w},p_{s}\\right)_{\\Omega} + \\langle\\mathbf{w}\\cdot\\mathbf{n}, p_{s}\\rangle_{\\Gamma} \\approx -\\left(\\nabla\\cdot\\mathbf{w},p_{s}\\right)_{\\Omega}\\; .\n",
    "\\end{align*}\n",
    "\n",
    "Complementary, the subscales mass conservation term is also integrated by parts as\n",
    "\\begin{equation*}\n",
    "\\left(q,\\rho\\nabla\\cdot\\mathbf{u}_{s}\\right)_{\\Omega} = -\\left(\\nabla q,\\rho\\mathbf{u}_{s}\\right)_{\\Omega} + \\langle\\left(\\rho q\\mathbf{n}\\right),\\cdot \\mathbf{u}_{s}\\rangle_{\\Gamma} \\approx - \\left(\\nabla q,\\rho\\mathbf{u}_{s}\\right)_{\\Omega}\\; ,\n",
    "\\end{equation*}\n",
    "as well as the subscales convective term\n",
    "\\begin{equation*}\n",
    "\\left(\\mathbf{w},\\rho\\mathbf{a}\\cdot\\nabla\\mathbf{u}_{s}\\right)_{\\Omega} = -\\left((\\rho(\\nabla\\cdot\\mathbf{a})\\mathbf{w} + \\rho\\mathbf{a}\\cdot\\nabla\\mathbf{w}),\\mathbf{u}_{s}\\right)_{\\Omega} + \\langle(\\mathbf{n}\\cdot(\\mathbf{a}\\otimes\\mathbf{w})), \\mathbf{u}_{s}\\rangle_{\\Gamma} \\approx -\\left((\\rho(\\nabla\\cdot\\mathbf{a})\\mathbf{w}+\\rho\\mathbf{a}\\cdot\\nabla\\mathbf{w}),\\mathbf{u}_{s}\\right)_{\\Omega} \\; .\n",
    "\\end{equation*}\n",
    "\n",
    "The subscale values are obtained from the finite element approximation residuals (Algebraic SubGrid Scales, ASGS) as\n",
    "\\begin{align*}\n",
    "\\mathbf{u_{s}} =& \\tau_{1}\\mathbf{R}^{M}(\\mathbf{u}_{h},p_{h})\\\\\n",
    "p_{s} =& \\tau_{2}\\mathbf{R}^{C}(\\mathbf{u}_{h},p_{h})\n",
    "\\end{align*}\n",
    "being $\\tau_{1}$ and $\\tau_{2}$ the stabilization coefficients\n",
    "\\begin{align*}\n",
    "\\tau_{1} =& \\left(\\frac{\\rho\\tau_{d}}{\\Delta t} + \\frac{c_{1}\\mu}{h^{2}} + \\frac{c_{2}\\rho||\\mathbf{a}||}{h} + \\frac{c_{3}}{\\kappa h^2}\\right)^{-1}\\; , \\\\\n",
    "\\tau_{2} =& \\mu + \\frac{c_{2}\\rho h||\\mathbf{a}||}{c_{1}} + \\frac{c_{3}}{\\kappa c_{1}}\\; .\n",
    "\\end{align*}\n",
    "being $c_{1}=4$, $c_{2}=2$, $c_{3}=2$ and $h$ the characteristic element size.\n",
    "Note that the selected subscales model contains higher order derivatives that vanish when using linear finite elements.\n",
    "Also note that the permeability $\\kappa$ can be alternatively expressed in terms of the resistance $\\sigma = 1 / \\kappa$.\n",
    "\n",
    "Thus, the final Galerkin residual weak form reads\n",
    "\\begin{align*}\n",
    "\\left(\\mathbf{w},\\rho\\mathbf{f}\\right)_{\\Omega} -& \\left(\\mathbf{w},\\rho\\frac{\\partial\\mathbf{u}_{h}}{\\partial t}\\right)_{\\Omega} - \\left(\\mathbf{w},\\rho\\mathbf{a}\\cdot\\nabla\\mathbf{u}_{h}\\right)_{\\Omega} + \\left(\\rho(\\nabla\\cdot\\mathbf{a})\\mathbf{w},\\mathbf{u}_{s}\\right)_{\\Omega} + \\left(\\rho\\mathbf{a}\\cdot\\nabla\\mathbf{w},\\mathbf{u}_{s}\\right)_{\\Omega} - \\left(\\nabla\\cdot\\mathbf{w},p_{h}\\right)_{\\Omega} - \\left(\\nabla\\cdot\\mathbf{w},p_{s}\\right)_{\\Omega} - \\left(\\nabla\\mathbf{w},\\mathbb{\\mathbb{C}}\\nabla^s\\mathbf{u}_{h}\\right)_{\\Omega} - \\left(\\nabla\\mathbf{w},\\mathbb{\\mathbb{C}}\\nabla^s\\mathbf{u}_{s}\\right)_{\\Omega}\\\\ \n",
    "&- \\left(\\mathbf{w},\\frac{\\mathbf{u}_{h}}{\\kappa}\\right)_{\\Omega} - \\left(\\mathbf{w},\\frac{\\mathbf{u}_{s}}{\\kappa}\\right)_{\\Omega} - \\left(q,\\frac{1}{c^{2}}\\frac{\\partial p_{h}}{\\partial t}\\right)_{\\Omega} - \\left(q,\\mathbf{a}\\cdot\\nabla\\rho\\right)_{\\Omega} - \\left(q,\\rho\\nabla\\cdot\\mathbf{u}_{h}\\right)_{\\Omega} + \\left(\\nabla q,\\rho\\mathbf{u}_{s}\\right)_{\\Omega} + \\langle\\mathbf{w}, \\left(\\mathbb{C}\\nabla^s\\mathbf{u}_{h}-p_{h}\\mathbf{I}\\right)\\cdot\\mathbf{n}\\rangle_{\\Gamma} = 0 \\; .\n",
    "\\end{align*}\n",
    "Note that the term $\\left(\\nabla\\mathbf{w},\\mathbb{\\mathbb{C}}\\nabla^s\\mathbf{u_{s}}\\right)_{\\Omega}$ vanishes if linear elements are used in the subscale approximation. \n",
    "\n",
    "## SYMBOLS TO BE EMPLOYED\n",
    "\n",
    "Shape functions $N_{I}$ and derivatives $\\nabla N_{I}$stored respectively\n",
    "in a vector $\\mathbf{N}$ and a matrix $\\mathbf{DN}$. Besides the\n",
    "following arrays are defined:\n",
    "\n",
    "- $\\mathbf{p}$ such that $p_{I}$ is the current step pressure of node I\n",
    "- $\\mathbf{p_{n}}$ such that $p_{n,I}$ is the previous step pressure of node I\n",
    "- $\\mathbf{p_{nn}}$ such that $p_{nn,I}$ is the two previous step pressure of node I\n",
    "- $\\mathbf{v}$ such that $\\mathbf{v}_{IK}$ is the current step velocity of node I, component K\n",
    "- $\\mathbf{v}_{n}$ such that $\\mathbf{v}_{n,IK}$ is the previous step velocity of node I, component K\n",
    "- $\\mathbf{v}_{nn}$ such that $\\mathbf{v}_{nn,IK}$ is the two previous steps velocity of node I, component K\n",
    "- $\\mathbf{w}$ such that $\\mathbf{w}_{I}$ is the linear momentum conservation test function value I at the considered Gauss point\n",
    "- $\\mathbf{q}$ such that $\\mathbf{q}_{I}$ is the mass conservation test function value I at the considered Gauss point\n",
    "\n",
    "Values on the gauss points are expressed in terms of such variables as\n",
    "\n",
    "- $\\mathbf{v\\_gauss}:=\\mathbf{v}^{T}\\mathbf{N}$ a 3x1 matrix\n",
    "- $\\mathbf{p\\_gauss}:=\\mathbf{p}^{T}\\mathbf{N}$ a 1x1 matrix\n",
    "- $\\mathbf{f\\_gauss}:=\\mathbf{f}^{T}\\mathbf{N}$ a 3x1 matrix\n",
    "- $\\mathbf{w\\_gauss}:=\\mathbf{w}^{T}\\mathbf{N}$ a 3x1 matrix\n",
    "- $\\mathbf{q\\_gauss}:=\\mathbf{q}^{T}\\mathbf{N}$ a 1x1 matrix\n",
    "- $\\mathbf{accel\\_gauss}:=\\left(BDF_{0}\\mathbf{V}+BDF_{1}\\mathbf{V}_{n}+BDF_{2}\\mathbf{V}_{nn}\\right)^{T}\\mathbf{N}$ a 3x1 matrix\n",
    "- $\\mathbf{pder\\_gauss}:=\\left(BDF_{0}p+BDF_{1}p_{n}+BDF_{2}p_{nn}\\right)\\mathbf{N}$ a 1x1 matrix\n",
    "- $\\mathbf{v_{h}:=v^{T}N}$ a 3x1 matrix\n",
    "- $\\mathbf{w_{h}:=w^{T}N}$ a 3x1 matrix\n",
    "- $p_{h}:=\\mathbf{p^{T}N}$ a 1x1 matrix\n",
    "- $q_{h}:=\\mathbf{q^{T}N}$ a 1x1 matrix\n",
    "- $\\mathbf{div\\_p}:=\\mathbf{\\nabla\\cdot p}$ a 1x1 matrix\n",
    "- $\\mathbf{div\\_v}:=\\mathbf{\\nabla\\cdot v}$ a 1x1 matrix\n",
    "- $\\mathbf{div\\_vconv}:=\\mathbf{\\nabla\\cdot vconv}$ a 1x1 matrix\n",
    "- $\\mathbf{grad\\_p}:=\\mathbf{DN^{T}\\cdot p}$ a 3x1 matrix\n",
    "- $\\mathbf{grad\\_v}:=\\mathbf{DN^{T}\\cdot v}$ a 3x3 matrix\n",
    "- $\\mathbf{grad\\_vconv}:=\\mathbf{DN^{T}\\cdot vconv}$ a 3x3 matrix\n",
    "- $\\mathbf{grad\\_q}:=\\mathbf{DN^{T}\\cdot q}$ a 3x1 matrix\n",
    "- $\\mathbf{grad\\_w}:=\\mathbf{DN^{T}\\cdot w}$ a 3x3 matrix\n",
    "- $\\mathbf{grad\\_sym\\_v}:=$ a 3x3 matrix (symmetric gradient in Voigt form)\n",
    "- $\\mathbf{grad\\_sym\\_w}:=$ a 3x3 matrix (symmetric gradient of the test function in Voigt form)\n",
    "- $\\mathbf{convective\\_term}:=\\mathbf{\\left(a_{h}\\cdot grad\\_vh\\right)^{T}}$ a 3x1 vector in where $\\mathbf{a_{h}}$ is the convective velocity. If the option ''Picard'' linearisation is set, the convective velocity is defined as a symbol $\\mathbf{v_{conv}}$ meaning that it is treated\n",
    "as a constant in the automatic differenctiation. On the contrary, if the option ``FullNR'' linearisation is set, the convective velocity\n",
    "is defined as $(\\mathbf{v_{h}}-\\mathbf{v_{mesh}})$ implying that $\\mathbf{v_{h}}$ is considered in the automatic differenctiation.\n",
    "- $\\mathbf{stress}:$ 6x1 shear stress vector computed by constitutive law using the previous iteration database. This stress vector is considered\n",
    "in the computation of the RHS, which is the residual of the previous iteration solution. When computing the LHS (RHS derivative w.r.t.\n",
    "the DOFs), the $\\mathbf{stress}$ is substituted by $\\mathbf{C*grad\\_sym\\_v}$ to consider the stress in the automatic differenctiation.\n",
    "\n",
    "## IMPLEMENTATION\n",
    "The residual functional implementation is splitted in two pieces\n",
    "- $\\mathbf{rv\\_galerkin}$: functional corresponding to the standard incompressible Navier-Stokes equations. If \"artificial\\_compressibility\" is set as True, the artificial compressibility terms are also added to the functional.\n",
    "- $\\mathbf{rv\\_stab}$: functional corresponding to the standard incompressible Navier-Stokes ASGS stabilization. If \"artificial\\_compressibility\" is set as True, the artificial compressibility stabilization terms are also added to the stabilization functional.\n",
    "- $\\mathbf{rv}=\\mathbf{rv\\_galerkin}+\\mathbf{rv\\_stab}$: summation of functionals. This is the functional used to compute the LHS and RHS elemental matrices.\n",
    "\n",
    "For the definition of the subscales, the residuals are splitted in\n",
    "- $\\mathbf{vel\\_residual}$: linear momentum conservation equation residual. \n",
    "- $\\mathbf{mas\\_residual}$: mass conservation equation residual. If \"artificial\\_compressibility\" is set as True, the artificial compressibility terms are also added to the residual.\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
